Evaporation and fluid dynamics of a sessile drop of capillary size 
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Theoretical description and numerical simulation of an evaporating sessile drop are developed. We 
jointly take into account the hydrodynamics of an evaporating sessile drop, effects of the thermal 
conduction in the drop and the diffusion of vapor in air. A shape of the rotationally symmetric drop 
is determined within the quasistationary approximation. Nonstationary effects in the diffusion of 
the vapor are also taken into account. Simulation results agree well with the data of evaporation rate 
measurements for the toluene drop. Marangoni forces associated with the temperature dependence 
of the surface tension, generate fluid convection in the sessile drop. Our results demonstrate several 
dynamical stages of the convection characterized by different number of vortices in the drop. During 
the early stage the street of vortices arises near a surface of the drop and induces a non-monotonic 
spatial distribution of the temperature over the drop surface. The initial number of near-surface 
vortices in the drop is controlled by the Marangoni cell size which is similar to that given by Pearson 
for flat fluid layers. This number quickly decreases with time, resulting in three bulk vortices in 
the intermediate stage. The vortices finally transform into the single convection vortex in the drop, 
existing during about 1/2 of the evaporation time. 



I. INTRODUCTION 



Evaporation of a liquid drop in an ambient gas is a well-known yet not completely solved problem of the classical 
physics (see, for example, P, Q)- Past decade was marked by significant advance in experiment and progress in 
understanding several key aspects of evaporation process, in particular, the vapor diffusion from the sessile drop 
surface and the hydrodynamic effects within the evaporating drops [1, |j, H, H, 0, H, Q ■ It was fomid, in particular, 
that the evaporating flux density is inhomogeneous along the surface and diverges on approach to the pinned contact 
line d, 01 . The resulting inhomogeneous mass flow modifies the temperature distribution over the drop surface and, 
hence, the Maraiigoni forces associated with the temperature dependent surface tension. The convection inside a 

droplet [1, i. El [nl m, [H, Q M, 

[1^ [TtI [TsI. [T9I appears quite different from the classical Marangoni convection 
in the systems with a simple flat geometry [2Cl |21|. Thermal conductivity of the substrate can also influence the 
formation of flows within a liquid drop since it is the magnitude of the conductivity which determines the sign of the 
tangential component of the temperature gradient at the surface close to the contact line and, therefore, the direction 
of the convection 7| . The observation of the distinct stages of the evaporation process [13, [H, [23| has revealed that the 
longest and dominating regime of the evaporation process is the constant contact area mode, where the contact line 
is pinned. As further the contact line gets depinned, the different regime, the constant contact angle mode switches 
on. Finally, the drying mode follows, in which the height, the contact area and the contact angle rapidly decrease 
with time. 

Understanding the details of the drop dynamics and evaporation is crucial to many ap plic ations involving this 
process, including preparing ultra-clean surfaces 25, 26, 27,[23|, protein crystallography [23, [30|, the studies of DNA 
stretching behavior and DNA mapping methods 31, 32, 33 1, developing methods for jet ink printing [H, [H, H^, and 



in many other fields (see, for example, [37|). The process of evaporation of the drop with the colloidal suspension in 
it is of interest for methods of fabrication of various structures on the substrate. One of the examples is the effect of 
evaporative contact line deposition, the so-called coffee-ring effect 0, 0, H, 38 , 39 , Ec | . Other important example is 
the self-assembly of long-range-ordered nanocrystal super lattice monolayer 4l|. |43. |43 1 . 

While the fundamentals of quasistationary evaporation behavior are established, to some extent [1, Q, the full 
dynamical description of the liquid drop evaporation is yet not available. In this work we undertake a step towards 
such a description offering a quantitative approach which enables us to account for all the relevant components of the 
evaporation process, namely, the fluid dynamics, the vapor diffusion, and the spatial temperature distribution in an 
evaporating sessile drop. The calculation are carried out according to the following scheme: (i) We derive the equation 
yielding the shape of the sessile drop and find how the evaporation rate is influenced by the deviations of the drop 
shape from the ideal spherical cap caused by the gravitational force, (ii) We solve the diffusion equation describing the 
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vapor kinetics and calculate the local evaporation flux from the surface of the drop as well as the resulting evolution of 
the drop evaporation rate with time. At this stage we take into account the nonstationary effects in a vapor diffusion 
and find that the dynamical corrections to the evaporation rate do not vanish exponentially, but decay as cx 1/ \/t. 
(iii) We solve the thermal conduction equation and find the spatial temperature distribution in the drop, (iv) We 
solve the Navier-Stokes equations and derive the velocity field corresponding to the convection within the drop, (v) 
Finally we arrive at a full description of the time evolution of the temperature and the fluid convection in the drop 
and identify characteristic stages of the thermocapillary convection. To this end we take into account inertial terms 
in the Navier-Stokes equations, including the time derivatives. 

To ensure self consistency of the derivation of the physical characteristics of the system we include the temperature 
variation over the drop surface (since the surface tension depends on temperature) into the boundary conditions for 
the fluid dynamics in the drop. Furthermore, we take into account that the velocity field influences the thermal 
conduction, since the velocities enter the thermal conduction equation. And, finally, we take into account that the 
local evaporation fiux is related to the heat transfer and, hence, to the temperature gradient on the drop surface. 

While the developed approach is quite general and thus applies to a wide variety of the evaporating problems, 
we focus here on the experimental situation corresponding to evaporation of the toluene sessile drop of the capillary 
size. The typical values of the characteristic parameters at Tq = 295K are taken from [i^l and are presented in the 
lower part of Table |T1 The initial values which we use in our calculations are given in the upper part of Table [J and 
correspond to the evaporating toluene drop containing the colloidal solution of gold nanoparticles. These nanoparticles 
stick to the drop surface making islands there and form a monolayer nanocrystal on the substrate after toluene dries 

out [Mm, Si- 

TABLE I: The parameter values used in the paper for modeling the evaporation and hydrodynamics of the drop. The tabular 
data are taken at To = 295K from 44]. 



Drop 
parameters 


Contact line radius 
Initial height 
Initial mass 


ro = 0.2 cm 
h = 0.1314 cm 
mo = 8.7 ■ 10"^ g 


Toluene 
characteristics 


Density 
Molar mass 
Thermal conductivity 
Thermal difFusivity 
Kinematic viscosity 
Dynamic viscosity 
Specific heat 

Thermal-expansion coefficient 
Surface tension 

Temperature derivative of surface tension 
Capillary constant 
Latent heat of evaporation 


p = 0.87 g/cm^ 

/i = 92.14 g/mole 

k = 1.311 • 10'^ W/(cm- K) 

K = k/{pcp) = 8.86 • 10"" cm^/s 

= 6.4 ■ 10"^ cmVs 
T] — up = 5.6 ■ 10^^ g/(cm- s) 
cv = 1.51 J/(g-K) 
P = 1.138 • 10"^ 
a = 28.3049 g/s^ 
a' = -da/dT = 0.1189 g/(s^- K) 
= 0.258 cm 
L = 300 J/g 


Toluene 
vapor 
characteristics 


Diffusion constant 

Saturated toluene vapor density 

Mean free path 


D = 0.1449 cmVs 

Us = 1.27 ■ 10-'' g/cm^ 

A = 2.7 • 10-^ cm 


Air 

characteristics 


Thermal conductivity 
Dynamical viscosity 
Density 

Kinematic viscosity 


ka = 0.258 • 10"^ W/(cm- K) 
r,a = 1.82 ■ 10-* g/(cm- s) 
Pa = 1.2 ■ IQ-^ g/cm^* 
I'a = r]a/Pa = 0.15 cm^/s 



We have identified three major distinct dynamic stages of the Marangoni convection during the evaporation of the 
toluene drop. During the initial stage, vortices appear near the surface of the drop, during the second stage they 
grow, coalesce and rapidly migrate into the drop bulk, and, finally, a single vortex survives governing all the fluid 
dynamics of the thermocapillary effect in the system. We derive the temperature profiles, evaporation rates, and the 
velocity distributions for each stage of the convection and uncover the important role of the convective heat transfer 
and the effect of inertial terms in the Navier-Stokes equations in the evaporation dynamics. 

The paper is organized as follows. In Sec. HI] we present basic equations and main analytical results along with 
the quantitative description of main dynamical processes. In particular, Sec. Ill Al deals with the shape of a sessile 
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drop. Various aspects of the vapor diffusion and the drop evaporation rates are discussed in Sec. Ill BHII Dl The 
characteristic Marangoni numbers for the drop are found in Sec. Ill El Sec lIIFl makes use of a self-consistent approach 
for calculating the velocity field and the spatial temperature. In Sec. IIIII we present our main results and discussion, 
including the experimental and theoretical data for the evaporation rates and contact angles, the dynamics of vortex 
structure and the temperature profile. Sec. IIVI contains discussion. The details of the developed numerical approach 
are described in Appendix El 



II. BASIC EQUATIONS AND FORMULAS 



A. The sessile drop shape 



During the evaporation process a drop loses its mass, hence with time its volume decreases and the shape changes. In 
this section we present the quasistationary method for calculating the shape of a sessile drop based on the hydrostatic 
Laplace equation. The approach holds as long as viscous forces, which generally enter the boundary condition for the 
pressure, are small. The ratio of viscous to capillary forces is characterized by the dimensionless number Ca = "rfv/a^ 
where v is the characteristic value of the velocity. In our case Ca « 2 • 10~* ^ 1. 

The Laplace equation states that the pressure difference Ap, taken at the different sides of the surface of the liquid 
in an arbitrary point equals crfc, where k — + l/i?2 is the mean curvature of the surface and a is the surface 

tension [i^ . Taking into account the gravity, one finds that at the surface of the drop 

k + = k ^ — ~ const, (1) 

where the shape of the drop surface is defined by the relation z = f{x, y) and a = y^2a/ (pg) is the capillary constant. 
The curvature of the surface of the drop is, in its turn, also expressed in terms of the function f(x, y). Indeed, 
k = —Tt{G~^Q), where G and Q are matrices of the first and the second quadratic forms of the surface [4y|. Plugging 
in the surface equation z = f{x,y) we find 

Q — I 1 + fx:fy \ Q ^ - I '^^^ '^^^ I (2) 

V fxfy ^ + fy J' i/i + + fl \ -f-y -fyy J 

Therefore, 



; rp '^fxfyfxy fxx{^ + fy) fyyi^ ~^ fx) 

We consider a sessile drop with the axial symmetry, where z — f{r), r ~ \/ + y^- Then Eq. ([3]) assumes the form 



(l + /'2)3/2 ^(l+J'2)l/2^ 



(4) 



where prime denotes the derivative with respect to r. It is convenient to parametrize and separate the variables 
as [47ll48t: 

r = r{s), z = z{s), (5) 

where s is the arc length on the drop surface taken from the apex. The curvature radii are Ri and i?2, where i?i is 
the curvature radius in the meridian rz-plane, and i?2 is the second curvature radius. Both curvature radii can be 
expressed in terms of the angle (j) between the normal vector to the drop surface and the symmetry axis. Namely, 

d(f> r 

^ = ^2 = ^— ■ (6) 

Hi els sm 

Also, dr/ds = cos0, dz/ds — — sincf) and r(0) = z(0) = (/)(0) = 0. Introducing the curvature radius Rq at the apex, 
we can rewrite the Laplace equation as 
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and therefore 



ds 



2 



pgz 

a 



sm ( 



(8) 



Hence, for the column vector y = {r{s), (j){s), z{s))'^ , where T is the transposition operation, we have the set of 
first-order differential equations with Cauchy boundary conditions: 



-T- = f(s,y), 
as 



f(s,y) = 

y(0) = (0,0,0) 



' Rq 

T 



P9± 
a 



f(o,y)= ( l,;^,0 



(9) 
(10) 

(11) 



where the last relationship ensures that (j){s)/r{s) ^ 1/Ro for s ^ 0. The Cauchy problem ([9|l- (fTT|) can be solved 
using the fourth order Runge-Kutta method or the Adams-Bashforth method [43|. Therefore, given the values of Rq 
and Smax, one finds the drop shape in the form of y(s) = (''(s), (t>{s), z{s))'^ . 
The quantities that are measured in the experiment are the mass of the drop 



m = TT/O 



r^(s) sin 0(s)ds , 



(12) 



and the radius rg of the substrate to which the drop is pinned. Given the values of Rq and rg, one can easily find 
Smax, and, therefore, m. Inverting the function m{Ro-,r) numerically, one can obtain the value of Ro and a drop 
shape for a given mass m and contact line radius tq. We note that Eqs. ([5|)- PT|) are exact for a sessile drop with axial 
symmetry. 

The evaporation rate for a given drop shape can be obtained as 



dm 




dt 


-L 



2TTr{s)J{r{s))ds, 



(13) 



where J(r(s)) is the mass evaporated per second from unit area of the surface, i.e. the local evaporation flux. 



B. Evaporation of the drop 

Consider a sessile droplet resting on a flat substrate. The vapor concentration above the droplet is time dependent 
and inhomogeneous during the evaporation process. Since the diffusion of the vapor from a near-surface layer is slower 
than the evaporation (45j , the vapor concentration at the surface of the drop is assumed to equal the saturation value 
Us- Far above the drop, the toluene vapor concentration is negligible. The dynamics of the vapor concentration in 
the surrounding atmosphere is described by the diffusion equation 



- = DAu. (14) 

We carry out our calculations (see Sec. IIIII and Appendix [A| taking into account both time dependence of vapor 
concentration and the deviations of the sessile drop shape from a spherical cap. It is instructive to begin discussing 
the problem with the more simple case allowing for an analytical description. To this end we notice that if evaporation 
can be viewed as an adiabatic process in a sense that the vapor concentration adjusts fast enough to the change of 
the drop size (and shape) i.e. on the time scales much less than the droplet evaporation time tf, then the diffusion 
equation can be replaced by the Laplace's equation Au = 0, and the evaporation kinetics can be considered as 
a steady state process. Indeed the time tp required for the Brownian particle to pass the characteristic length vq is 
tp = Tq/D ^ tf under the experimental conditions tp « 0.2 s and t/ ~ 500 s. At the same time it interesting as well 
to consider particular dependencies of the local evaporation rates on time. Our simulation results based on Eq. (|14|) 
for the drop with the fixed surface can be fitted with the power-law time dependence 

J(M)^J(.c„)(l + ^), (15) 
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which is almost exact with the accuracy within 1% for t > 0.5 s and |ro— r| > 0.01 cm (see Fig.[T]). Here J(r, t) = |DVu| 
is the local evaporation flux on the drop surface, and constant A = 0.966 is the only fitting parameter. For the case of 
the spherical cap the asymptotic value of the evaporation flux density J(r, oo) can be related to J(^, 9) from Eq. ([20|) . 
Eq. (|15p and Fig. [T] confirm the above estimation that the vapor concentration becomes stationary under the condition 
t ^ rg/D. At the same time, ioi t = tf the second term in (fTS]) still exceeds 1% of the first term. The time dependence 
J(r, t) — J{r, oo) cx in Eq. (jlSP is known for an isotropic diffusion when the vapor concentration is kept saturated 
on the fixed surface of the sphere Q . This is also valid for a diffusion from a flat plane [45^ . Our results demonstrate 
that such dependence also takes place for an inhomogeneous diffusion from a fixed surface of the sessile drop. 

We further derived the numerical results which take into account the time dependent drop profile (see Fig. [5]). In 
particular, they show that the corrections for the evaporation rate due to the nonstationary effects may be up to 5% 
of the resulting value. 



FIG. 1: (Color online.) Local evaporation flux density J{r,t) from a fixed surface as a function on r for t = 
1,2,3,5,10,50, 100,500, 1000,4000 seconds of the evaporation process, from top to bottom, respectively (black curves). The 
time dependence of the vapor concentration is taken into account, whereas the surface of the drop with m = 8.7 mg is fixed. 
The surface is taken as a spherical cap (left panel) and as a sessile drop surface (right panel). The bottom line on the left panel 
(blue curve) represents J(r, oo) taken from the exact solution ([20}. The bottom line on the right panel (green curve) represents 
J(r, co) for a sessile drop. 




005 0.1 0.15 0.2 ' 0.05 0.1 0.15 0.2 ' 



Deegan et al.[3] reported an analytical solution for a stationary spatial distribution of the vapor concentration for 
a drop with the shape of a spherical cap. The problem was shown to be mathematically equivalent to that solved 
by Lebedev [i^, who obtained the electrostatic potential of a charged conductor having the shape defined by two 
intersecting spheres. Modifying the equations derived in [J |4^, we arrive at the following formula for the vapor 
concentration u(a, (3) above the surface (pS)) of a spherical drop: 

, . , , I— ; f°° coshdr cosh/3r 

u(a,/J) = Uoo + (us - Uoo)v2(coshQ; - cos/3) / —. -— P_i/2+ir (cosh Q;)dT. (16 

Jq cosh ttt cosh(TT — 0)t ' 

Here 6 is the drop contact angle, — tt + 9<f3<Tr — 9, Uoo — 0, a and (3 are the toroidal coordinates: 

ro sinh a tq sin f3 



cosh a — cos (3 ' cosh a — cos P 

At the drop surface P ^ n — 9, therefore £, — t/tq = Va;^ — 1/ {x + cos 9), where x(^, 9) = cosh a, hence, 



The evaporation flux at the surface is 



D du 
J — IDVmI — — (cosh a — cos/3) —— 
ro 9/3 



and, therefore. 



(17) 



x{^,9) ^cosha = ^—^ . (18) 



(19) 



Ji^,0)^ — ('^ + V2(x(e, 9) + cos 9)'/' r tanh(7r - 9)tP_,/,+„ (x(C, 9))dr\ . (20) 

ro V 2 Jo cosh TTT J 
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After substituting ^ = r/rp and (fTSl) in Eq. ((20)) . one obtains J(»'). Fig. [2] shows log J(r) for drops with different 
values of the contact angle 9. 



Log J(r/R) 
5.75 h 




|Log(1-r/R)| 



FIG. 2: Log-log plot of J(r) for drops with different contact angles. Lines marked with a, b, c, d, e correspond to values vr/f 20, 
it/6, n/'i, 7r/3, 237r/60 of the contact angle 6, respectively. 



The total drop mass is 



and the total evaporation rate is 



purl 



„ tan - 3 
6 2 



tan^ - 



dm 
'dt 



27rr5 



- sin^ 



(21) 



(22) 



While it were possible to simplify the expression for the drop evaporation rate adapting the sessile drop shape to 
that of the spherical cap of the same mass (see Sec. Ill D[) . and then integrate Eqs. ipDl) . ([H]), we will employ a more 
general method and develop a description of the evaporation rate which takes into account altogether spatial and time 
variations of the vapor diffusion and the sessile drop shape. 



C. Estimation for the evaporation time 



In this section we derive estimates for the evaporation time of the drop, which in spite of their relative simplicity, 
offer, nevertheless, a fairly good description of the experimental results. Let us consider a spherical evaporating drop 
of radius ro = 0.2 cm, so that the vapor density is saturated at the drop surface and vanishes far away from the 
drop. The diffusion of the toluene vapor controls the evaporation process. The outward flow of the vapor through 
the concentric sphere around the drop of the radius R is J = —D ■ AirR^du/dR and does not depend on R. Here u 
is the vapor density on the surface of the sphere of radius R, and D is the diffusion constant for the toluene vapor. 
Therefore, R^du/dR = —A, where A = J/{ATrD). Since the vapor concentration far away from the drop is negligible, 
we obtain u = A/R, A = vqUs and J = ATrDroUg. On the other hand, the flow of molecules from the fluid drop is 
J = —pd{ATrrQ/3)/dt. Therefore, rgdro/dt = —Dus/p, and the evaporation time of the drop can be written as 

t^JLlL (23) 
Us 2D' ^ ' 

Eq. is a well-known result of the classical Maxwell's theory of the evaporation (see, for example, [2]). The diffusion 
constant for toluene vapor can be estimated as 

D^^{v)X. (24) 



Here (v) = y/8RT/ (ttp) = 26035 cm/s is the average velocity for the thermal motion of the vapor molecules, A is the 
mean free path which may be estimated as follows. We replace the toluene molecules by the impermeable balls of the 
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diameter d ~ 1 nm (we found the estimate for the diameter by adding together the lengths of the chemical bonds and 
taking into account the molecule geometry). Therefore, the mean free path is A = \|{'K^/2n£■) = kT/{TTV2Po(f) « 
2.7- 10~^ cm. Here n is the density of the toluene molecules in the toluene vapor. If we substitute the obtained values 
into the Eq. ([M]) . we get the estimate D ~ 0.23 cm^/s. This qualitative estimate is about in 1.5 times larger than the 
result obtained from comparison of our numerical results and the experimental data in Sec IIIII 

Substituting the value of D to the Eq. we get the evaporation time of the spherical drop: t « 594 s. 

D. Approximating a sessile drop surface with spherical caps 

Spherical caps offer a very simple model allowing for carrying out analytical calculations to the end and are thus 
widely used for the modeling of the droplet. The equation for the surface is 

Figure [3] shows surface of the sessile drop and three approximating spherical caps with the parameters listed in 
Table HTD] The radius of the substrate is tq — 0.2 cm for all the drops. The first spherical surface has the contact 
angle of the sessile drop. The second spherical surface has the mass of the sessile drop. The third spherical drop has 
the height of the sessile drop. 

As seen from the Fig. |3] and Table WTPj the characteristics m, h and 6 of all three approximating spherical caps are 
quite close to those of the sessile drop. This is not the case for local parameters such as the local curvature of the 
sessile drop. As follows from Eq. ([1]), for the local curvature to be approximately constant along the drop surface, the 
condition Bo = pghro/ {2a sin 6) <C 1 has to be satisfied. Here Bo is dimensionless number which is analogous to the 
Bond number. 

FIG. 3: (Color online.) The sessile drop surface (red curve) and three spherical caps (curves 1,2,3; see Table HT^ . 



z 




TABLE II: Characteristics of the sessile drop and three spherical caps. 



Parameters 


sessile drop 


cap 1 


cap 2 


cap 3 


m, mg 


8.7 


9.95625 


8.7 


8.22035 


6, radian 


1.3032326 


1.3032326 


1.20453564 


1.16292249 


h, cm 


0.13145179 


0.152552 


0.137494 


0.13145179 


X{e) = (7r-26l)/(27r-26l) 


0.145545 


0.145545 


0.189081 


0.206135 


curvature at top, cm~^ 


8.0611 


9.64418 


9.33673 


9.17966 


curvature at the contact line, cm~^ 


12.0727 


9.64418 


9.33673 


9.17966 



For our sessile drop Bo ~ 0.4; thus its curvature changes over the surface by factor of 1.5. Furthermore, one sees 
from Fig. 3] and Table IIIDI that the assumption of the spherical cap shape cannot provide an approximation for 
the local evaporation flux density J(r), which would have been accurate enough over a whole surface. At the same 
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FIG. 4: (Color online.) Black curve - J(r) for the sessile drop with m — 8.7 mg. Red, blue and green curves - J(r) for the 
first, the second and the third spherical drop corresponding to notations of Table HT^ 



time, it is not surprising that the total evaporation rate, integrated over the surface of the drop, does not depend 
much on whether the drop is sessile or spherical. With the values from Table HI one may numerically integrate ((22|) 
and obtain the evaporation rates for the three spherical drops: dmi/dt = 20.52 /^g/s, dm2/dt = 19.76 /ig/s, 
dm^/dt — 19.47 /ig/s. The difference is within a few per cent. 



E. The Marangoni convection 



A fluid convection within a drop caused by the temperature dependent surface tension which persists in the drop 
during the all stages of an evaporation process is called the Marangoni convection, and was first observed in its classical 
form by Benard in a process of formation of the characteristic hexagonal convection patterns in flat fluid films. A 
theoretical description of the Marangoni convection was developed by Pearson [2l| . The Marangoni number Ma (see 
Eq. (j26p below) characterizes a relative importance of surface tension forces caused by the temperature variation and 
viscous forces, leading to the Benard-Marangoni instability if Ma exceeds a critical value M^. For a flat fluid film, 
Mc « 83 (2l| . To estimate the Marangoni number for a drop we exploit a formal similarity between a drop and a 
flat layer. We extract the temperature difference between the substrate plane and the apex of the drop from our 
simulation results as AT « 1 K. The drop height h, the thermal conductivity k, the kinematic viscosity ly, the specific 
heat cv, and the rate of the change of the surface tension with the temperature, —da/dT, are taken from the TablelH 
Then the value of the Marangoni number at the beginning of the drop evaporation process is 



Ma = "^^^ \ — w 2800. 26 

i^k 

According to experimental data, a turbulent regime arises for considerably lar ger values of the Marangoni number. 
For example, in the water drop the turbulence takes place for Ma > 22000 [50|. As the Reynolds number for 
toluene drop. Re « 62, is also comparatively small, the flow in our problem is of a laminar character. A competition 
between the buoyancy-induced flow and the thermocapillary convection is characterized by the dimensionless number 
B — pgh'^ [3 / (J .14,a') [211. In our case B « 0.02 ^ 1, and the buoyancy- driven convection is negligibly small. 

For a flat liquid layer the height of the Marangoni cell coincides with the layer thickness /i, whereas the transverse 
size of the cell is determined by the characteristic instability wavelength A. If the Marangoni number well exceeds the 
critical value, then A = 27r/i/a and a = ^/Ma/S pll |. 
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Our simulations show that near the surface a "street" of vortices forms at the initial stage of the convection (see 
Sec, imp . The size of the vortices (i.e. the distance over which they extend into the bulk of the drop), hd, constitutes 
a noticeable fraction of the height h of the drop. With the assumption that the instability wavelength A and the 
Marangoni number can be described similarly to the case of flat liquid layer, one can take the value hd for the layer 
thickness and obtain A = 2Tr\/hhd/a. Taking as an example the number of the near-surface vortices A'^ = 9 and 
Smax = 0.256 cm, we find A = s^ax/N sa 0.028 cm and 

This although crude estimate turns out to be in a good agreement with our simulations for the Marangoni convection 
(see Figs. 1617151 and their discussions). 

F. Hydrodynamic equations 

The basic equations inside the drop are the Navier-Stokes equations, the continuity equation for the incompressible 
fluid, and the equation for thermal conduction 

dv 1 

— + (vV)v + - gradp ly Av, (28) 
divv = 0, (29) 

dT 

— +vVT = kAT. (30) 

Here A = d"^ jdr^ + d/rdr + d'^ /dz^, v = r\l p is kinematic viscosity, k — k/{pcp) is thermal diffusivity. Other terms 
in the thermal conduction equation i^/(2cp) {dvi/dxk + dvk/dxi)'^ , where sum over i and k is assumed, are negligibly 
small. 

Applying the "rot " operation to the both sides of equation (|28p , one excludes the pressure p and obtains 

d 

—rot V + (vV)rot v — (rot v • V)v — i/Arot v. (31) 



Therefore, 



dt 

where the vorticity 7 is introduced by 



|-7(r, z) + (vV)7(r, z) ^ ( Aj{r, z) - ) , (32) 



7(^2) = ^-^; rotv = 7(r,z)iy. (33) 
oz or 

We define the stream function i/), such that 

dip dip 

— =rvr : -TT = -i"Vz (34) 
oz or 

One has d'^ip/{drdz) = Vr + r{dvr/dr) = —r{dvz/dz), therefore Vr/r + dvr/dr + dvz/dz — 0, i.e. the velocities 
obtained with Eqs. fM]) will automatically satisfy the requirement 

Applying the Laplace operator to the stream function, one obtains Aijj — — 2vz , which can be transformed to a 
more convenient form Aij) = r"/ with the modified operator A that differs from the Laplace operator by the sign of 
the term d/{rdr): 

AiP=^-l^ + ^^r(—-—']=r'y (35) 
dr^ r dr dz'^ \ dz dr J 

The method of the numerical solution is as follows: at each step, inside the drop 

1. We solve Eq. with the proper boundary conditions and find 7(r, z); 
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2. We solve Eq. p5|) and find -(/'(r, z). Then we find velocities with Eqs. p4)) : 

3. We solve Eq. (l30|) and obtain the new boundary conditions for Eq. ((32|) . 

Details of the numerical approach are given in Appendix [X] The proper boundary conditions for quantities 7 and 
satisfying Eqs. and (P5|) respectively, are derived in Appendix [Bl They take the form 7 = 0forr = 0;7 = dv^/dz 
for z = 0; 7 = da /{rids) + 2vrd(j)/ds on the surface of the drop; ip = at all boundaries: at the surface of the drop, 
at the axis of symmetry of the drop (r = 0), and at the bottom of the drop (z = 0). 

Here da/ds = —a'dT/ds is the derivative of surface tension along the surface of the drop, where the distribution 
of temperatures, which is found with Eq. ([50)) is taken into account, and 

a{T)^ao-a'{T~To), (36) 

is the experimental dependence of surface tension a on temperature T. 

The boundary conditions for the Eq. §UI are dT/dr = for r = 0; T = To for z = 0; 

dT/dn = -Qo(0/fc (37) 

on the surface of the drop. Here Qa{r) = LJ{r) is the rate of heat loss per unit area of the upper free surface, J(r) is 
the local evaporation flux, Tq is the temperature of the substrate, and n is a normal vector to the surface of the drop. 

The relation Q{){r) — LJ{r) implies that the heat flow from ambient air towards the drop surface is negligible. 
This is the case if the temperature difference between the drop surface and the air far from the drop is less than 
This is well satisfied in the problem in question. 

III. RESULTS 
A. Evaporation rates 

The evaporation rate of sessile drops is controlled mainly by the vapor diffusion in the surrounding atmosphere . 
Here we determined the diffusion coefficient of toluene vapor in air by measuring the toluene evaporation rate. The 
main conditions and parameters of the experiment are as follows. The sessile drop of lOfxl of toluene is lying on a 
substrate and evaporates to the atmosphere during « 500 s. The contact line of the drop is pinned to the edge of the 
substrate. The weight of the drop was measured every 10 seconds. As substrates, we used silicon wafers with a 100 
nm thick amorphous silicon nitride layer. 

The results of the measurements of drop evaporation rates are presented in FigO Open circles are experimental 
values for the pure toluene evaporation and open triangles are values for the gold nanoparticles colloid. The solid line 
is the simulation result obtained in terms of Eq. (|13p within the full numerical scheme. As seen in Fig. [51 it agrees 
well with the experimental data. Comparing the experimental data and the simulation results permits us to find the 
diffusion coefficient of the vapor, D, the only parameter controlling evaporation, a.s D = 0.1449 cm^/s. It is worth 
noting that the rough estimate mentioned in Sec. Ill Cl is larger by the factor of about 1.5. 

One can divide the evaporation process into three characteristic parts. Our simulation describes the main part, 
i.e. the quasi-stationary diffusion-limited evaporation, which lasts for about 400 seconds. While there is a good 
agreement of the simulation and experimental data during the diffusion-limited regime of evaporation, this is not the 
case for the drying stage, i.e. the final 50-100 seconds. At the drying stage, the experimental evaporation rate drops 
rapidly but not abruptly. The main physical reason for such a behavior is the depinning of the contact line. The 
interface area then shrinks over time. The two distinct parts in the drying stage of evaporation can be identified. 
The first part, that occurs in the i=400— 510 time window can be described as din/dt(x{to—t)" that fits well the 
experimental data with io=550(2), a=0.41(2) for the evaporation of the pure toluene and ^0=522(3), a=0.28(8) for 
the colloid evaporation. The last part of the drying stage is a pronounced exponential decay of the evaporation rate 
dm/dtocexp^—t/td), with the decay time td— 20.5(5) for the colloid evaporation. The exponential behavior during 
this stage qualitatively follows from the fact that change of the drop surface area A controls here the evaporation rate 
evolution: dA/dt oc dm/dt cx A. Inset in Fig. [5) represents the experimental mass variation during the drying stage of 
the evaporation process. 

B. Street of vortices and the dynamics of vortex structure 

Our calculations demonstrate the presence of the characteristic early regime in the dynamics of the Marangoni 
convection in an evaporating sessile drop. For various liquids and drop sizes, the vortices arise near the surface of the 
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FIG. 5: Drop evaporation rate dm/dt. Experimental data for the pure toluene evaporation (open circles) and An nanoparticle 
colloid (open triangles). The solid line describes simulation results obtained within the full numerical scheme for the vapor 
diffusion (see also Sec. IIIB|I . The dotted line is obtained assuming the local evaporation flux to be uniform over the drop 
surface and constant with time. Inset: Mass variation at the end of evaporation. 
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drop. For a toluene drop, this regime quickly arises and evolves up to i ~ 0.3 s. This is quite a short time period 
as compared with the total evaporation time ~ 550 s, but it admits an experimental study. The vortices grow, the 
number of vortices decreases, and eventually they evolve into the Marangoni convection cells in the bulk of the drop. 
This can be seen in Figs. I6|71 where the vortex structure contains six pairs of near-surface vortices and a corner vortex, 
and the temperature displays just six humps at the surface. As was shown in Sec. Ill El this number of vortices in 
the drop is controlled by the Marangoni cell size which is similar to that given by Pearson for flat fluid layers. The 
existence of near-surface vortices and the associated humps in the profile of the surface temperature, become more 
pronounced with the decrease of the viscosity of the liquid. There are no near-surface vortices when the viscosity 
increases more than in four times as compared with the toluene drop. 

As seen in Fig. [71 the extrema of the surface temperature correspond to the change of sign of the tangential 
component of velocities at the surface. The reason for this behavior is that the fluid flow moves from the higher to the 
lower temperature regions of the surface, because, according to (|36p . the surface tension decreases with increasing the 
temperature. The flows result in a redistribution of the temperature due to the convective heat transfer. The number 
of the surface temperature humps and the number of the near-surface vortices decrease during their evolution. 

The initial conditions chosen when generating Fig.[6]were: the vanishing velocities, constant room temperature, and 
the vanishing vapor concentration. For a description of a real experiment they have to be slightly modified to include 
the weak stochastic distribution of the surface temperature and velocities. We have carried out such calculations 
using a small-scale stochastic initial conditions. This has changed the particular behavior of the fiuid dynamics 
only on the initial stage of the process, where a large number of small-scale surface vortices arise. Then the surface 
vortex structure quickly evolves into exactly the same one as we obtained for the basic initial conditions. This result 
demonstrates the generic character of the near-surface vortex regime in Fig. [6] at the early stage of the formation of 
the Marangoni convection. 

An initial velocity field within the sessile drop can also be strongly disturbed right after the drop has fallen down 
on the substrate, or for some other reason. We model such a situation by choosing random initial conditions for the 
bulk velocity field, which are in agreement with the continuity equation (j29[) . We find that strong disturbances of the 
bulk velocities can noticeably modify the initial stage of the drop dynamics, but they do not modify its main stage. 
For the initial random velocities of the order of 5 cm/s (which well exceeds the typical velocity in the vortex 1 cm/s), 
the difference between the dynamics of the disturbed and resting drops disappears after t « 0.5 s. This verifies the 
stability of the large-scale drop dynamics with respect to disturbances of the initial temperature and velocity fields. 
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FIG. 6: (Color online.) The velocity field at t = 0.16 s. 





During the enlargement of the near-surface vortices, their number decreases, and the convection involves the bulk 
of the drop. As a result, for t « 0.45 s, three bulk vortices control the velocity and temperature fields in the drop, 
as seen in left panels in Figs. I8I9I During the coexistence of three vortices, the corner vortex starts growing at the 
expense of the other two vortices, and eventually at t « 2.0 s it occupies the whole drop volume. A spatial dependence 
of the temperature along the drop surface is nonmonotonic, if the drop contains more than one vortex (see Fig. [9]). 
Right panels in Figs. 18191 demonstrate how in the single- vortex regime effects of Marangoni forces drive liquid along 
the surface to the apex, where the fluid penetrates along the symmetry axis in the depth of the drop. 
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FIG. 8: (Color online.) Left panel: The velocity distribution at t — 0.5 s. The stage of drop dynamics with three vortices takes 
place from t ~ 0.45 s to t ~ 2.0 s. Right panel: The velocity distribution at t = 30 s. A distribution with a single vortex takes 
place from t « 2.0 s to t ~ 250 s. 
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FIG. 9: (Color online.) Distributions of temperature within drop for the stages of drop dynamics with three vortices (left) and 
a single vortex (right). The distributions are taken eXt = 0.5 s and at f = 30 s. correspondingly. Temperature scale shown at 
the right column. 




The regime with the single vortex represents one of the main stages of the dynamics of the evaporating sessile drop. 
It lasts up to t w 250 s. More than half of the drop mass evaporates during this period of time. If the initial values of 
mass, height and contact angle of the drop are m = 8.7 mg, h — 0.1314 cm, B — 1.2045, then at the moment t = 250 s 
we find m = 4.0 mg, h = 0.0685 cm, Q — 0.716. In particular, h/{2rQ) w 0.17, i.e. the drop shape is noticeably 
flattened. The total time of the evaporation is 508 s. 

The quasistationary single-vortex state loses its stability dX t k, 250 s and the vortex acquires a pronounced non- 
stationary character. During this nonstationary regime, the fluid pulsations take place. The characteristic frequency 
of the pulsations corresponds to the circulation period 0.15 s of a fluid element in the original vortex. Initially the 
pulsations are concentrated near the center of the original vortex. Then, as shown in Fig. 1101 the single-cell pulsating 
state breaks into two-center (and later three-center) pulsating structure. Eventually at t ~ 300 s a quasistationary 
state with three vortices arises. 

The numerical calculations of the fluid dynamics were tested with several different mesh sizes. The respective 
results are qualitatively identical and show reliable convergence of the quantitative characteristics. For example, the 
single-vortex regime was found to arise at 3.48 s, 2.48 s, 2.2 s, 2.06 s, 2.01 s for 100x100, 150x150, 200x200, 250x250, 
300x300 mesh elements covering a half of the drop cross-section. 
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FIG. 10: (Color online.) Fluid pulsations ai t — 263 s. (left panel) and at t = 281 s. (right panel). 




FIG. 11: Temperature profile along the surface as function of r, in disregarding the effect of fluid flow on the thermal conduction. 




C. Temperature profile 



If the thermal conductivity of a substrate is large compared to that of the liquid, then the temperature can 
be maintained practically constant at the substrate-fluid interface. This is the case, in particular, for the silicon 
nitride substrates used in experiments (4ll . 14^ . |43| . The silicon nitride is a material with high thermal conductivity, 
approximately in three orders larger than for the toluene. For this reason, the boundary condition for the temperature 
distribution at the substrate can be reduced to the constant temperature. Heat transfer between the substrate and 
the drop plays an important role in establishing the temperature profile in the drop. This also excludes a possibility 
for the reversal of the Marangoni convection [f\ , taking place for substrates with relatively small thermal conductivity. 

The characteristic scale of the temperature variation on the drop surface can be easily estimated if one disregards 
the effect of velocities on the thermal conduction. According to the evaporation rate data, the heat loss per unit area of 
the drop near apex is Qq = LJ{r = 0) w 0.027 W/cm^ during the early stage of the evaporation process. Therefore, 
as follows from the boundary condition (pT]). we obtain \dT/dn\ = 22.88 K/cm. The calculations show that the 
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temperature dependence is almost linear on z at r = in disregarding the velocity term in thermal conduction. This 
permits to obtain the temperature difference between the substrate and the apex of the drop as 6T « \9T / dn\^^Q ho = 
S.OOSTi^. This estimate is in a good agreement with the temperature profile along the surface shown in Fig. [Til where 
the effect of the velocity field on the temperature in Eq. ([50)1 is disregarded. The temperature variation in Fig. [TT] has 
larger amplitude and takes the monotonic form, as compared with the respective curve in Fig. [T] This demonstrates 
that the fluid flow noticeably modifies the temperature variations in the drop. 

Similar consideration for the evaporating water drop studied in results in Qq = 0.093 W/cm^, \dT/dn\ = 
15.49 K/cm, ST « 0.56 K. This estimation for the temperature difference exceeds in about 24 times the respective 
numerical results obtained in Q (see Figs. 1 and 2 there). We consider this discrepancy as an internal contradiction 
of Q leading to a significant underestimation of the velocity values in the water drop. Our numerical calculations 
confirm that ST « 0.6 K for the water drop of At the same time, our calculations exactly reproduce the velocity 
field shown in Figs. 4 and 5 in Q, if one takes a given surface temperature distribution of Q. One should also note 
that the role of convective term in thermal conduction equation is substantially less for the water drop of due to 
a small size of the drop. 



D. Contact angles 



The evaporation process with the pinned contact line is accompanied with an increase of the oblateness of the drop 
shape. Figure [T^] displays our results for time dependence of the contact angle of sessile drop. Our results agree 
with the time-dependent contact angle determined experimentally under the identical initial drop parameters for the 
colloidal solution of gold nanoparticles in the toluene drop [42] . In [42] highly ordered nanoparticle islands were found 
to form on the drop surface. For this reason their orientation angle, which was determined experimentally close to 
the contact line with small angle X-ray scattering, coincides with the drop contact angle. 




IV. DISCUSSION 



We have developed the approach for studying the evaporation and fiuid dynamics of a sessile drop of a capillary 
size and applied it for the description of the toluene drop evaporation. In |.41i . ,42. . ,43j the evaporating toluene drop 
containing colloidal solution of gold nanoparticles was used for realizing the self-assembly of nanocrystal monolayer. 
During the evaporation process, the nanoparticle islands were found to form on the drop surface [42l ]. For under- 
standing the island formation as well as nanoparticle segregation, a study of the fiuid and evaporation dynamics is 
necessary. In particular, the convective flows described in this paper are important for driving the nanoparticles to 
the contact line and the drop surface. 

The numerical simulations were carried out and the system of the diffusion equation for the vapor, the thermal 
conduction equation for the temperature and the Navier-Stokes equations for the fiuid flow in the drop is solved. The 
shape of the drop is controlled by the quasistationary Laplace equation, which includes the effect of the gravitational 
forces. Their role in forming a proflle of a nonspherical sessile drop is characterized by dimensionless number Bq = 
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pghro/ {2a sin 9), which is analogous to the Bond number. Since in our case Bo « 0.4, effect of gravitational forces, 
in general, should be taken into account. We have found that deviations of the sessile drop shape from the spherical 
cap are noticeable in local drop characteristics such as the local curvature and the evaporation flux density. At the 
same time the integrated over the drop surface characteristics, like the rate of the mass loss, are well described by the 
spherical cap approximation. 

The experimental and simulation results for the time dependent drop evaporation rate agree well during the main 
longest stage of the evaporation process. The calculated evolution of the contact angle agrees with the results of 42]. 
In studying nonstationary diffusion equation, we found that the time-dependent corrections to the stationary local 
evaporation rates of the sessile drop do not vanish exponentially, but decay much slowly cx 1/ Vt with time and, 
in general, can be noticeable. This kind of behavior was established earlier only for homogeneous diffusion from a 
surface. 

Our solution demonstrates the presence of several time stages in the evolution of the Marangoni convection in 
the drop. The main quasistationary fluid flow, containing only one vortex in the drop, is formed during the longest 
period of the evaporation process. Disregarding the time derivatives and the dynamical behavior of the quantities 
would result directly in a stationary state with one vortex in the sessile drop for a given contact angle. This kind of 
single- vortex solutions was obtained in previous theoretical studies of the sessile drop evaporation [a . Taking into 
account an explicit time dependence of all quantities, permitted us to identify a sequence of early dynamical stages of 
the Marangoni convection in the drop, containing a street of near-surface vortices which then transforms to the state 
with three bulk vortices. The stability of the large-scale drop dynamics with respect to perturbations of the initial 
conditions is identified. Also, for a flattened evaporating drop we found under the conditions in question a three- vortex 
state instead of the single-vortex one. We establish an important role of inertial terms in the Navier-Stokes equations 
and convective heat transfer terms in the thermal conduction equation. 

Our further plans in developing the present approach include investigations of effects of surface surfactants and 
nanoparticles, effects of a nonstationary dynamics of the drop profile and a detailed study of the infiuence of substrate 
properties on the drop dynamics. 

Authors thank H.M. Jaeger for fruitful and encouraging discussions. This research was supported by the U.S. 
Department of Energy Office of Science under the contract DE-AC02-06CH11357 and by the Program of Russian 
Academy of Sciences. L.Yu. Barash thanks the Dynasty Foundation for the partial support of his research. 

APPENDIX A: NUMERICAL METHODS 
1. A brief outline of the method 

The simultaneous calculation of the physical quantities in the drop can be partitioned into several steps: 

1. We apply to the diffusion equation (jl4|l the implicit finite difference method using irregular mesh outside the 
drop and a variable time step. We use a boundary interpolation in a vicinity of the drop surface. For the 
boundary conditions we take u = Us on the drop surface, u = far away from the drop, du/dr = and 
du/dz = on the axes r = and z = correspondingly. 

2. Calculations of the stream function ■;/; and velocities v inside the drop are based on Eqs. p4|) . ((35)) . The implicit 
finite difference method with a regular mesh inside the drop is applied. We use a boundary interpolation near 
the drop surface. For the boundary conditions we take '0 = at all boundaries: on the surface of the drop and 
on the axes r = and z = 0. 

3. We solve Eq. (|32p to obtain the vorticity 7 inside the drop. The explicit finite difference method with a regular 
mesh inside the drop is used. We use a boundary interpolation close to the drop surface. For the boundary 
conditions we take 7 = for r ^ Q\ ^ = dv^jdz for z = 0; 7 = da/ (j]ds) + 2vTd(j)/ds on the drop surface, where 
da/ds — —a'dT/ds is the derivative of the surface tension along the drop surface according to psp . 

4. For calculating the temperature T inside the drop, the explicit finite difference method with a regular mesh is 
applied to the thermal conduction equation (|30p. We use a boundary interpolation in a vicinity of the drop 
surface. The boundary conditions take the form dT/dr = for r = 0; T = To for z = 0; dT/dn = —Qo{r)/k = 
—LJ{r)/k on the drop surface. Here Qo{r) is the rate of heat loss per unit area of the free surface, n is a normal 
vector to the drop surface. 

5. During the iterative procedure, the drop shape is recalculated in accordance with the evaporative mass loss for 
the respective time interval. For this purpose we use the Runge-Kutta method for Eas.([9))- pT|) . 
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2. Explicit method for temperatures inside the drop 

For solving Eq. ((30|) the following explicit finite difference method is used: 



rpn-\-l rpn rjin T^n rpn rpn 



which gives 

T^+' = a.{Tl^+, + 77:,-_i) + a,.(l + l/(2z))7^"+i^, + a,.(l - l/(2*))7^"_i,, + (1 - 2a, - 2a,)T^- 

— PrVrij{T^+lj — 2T-l,i) ~ /^zf zij — T"j^i)- (A2) 

Here T^" is the temperature at the n-th time step, i,j are the coordinates on the regular mesh, hr,hz are steps along 
the respective axis, ht is the time step, 

ar^Kht/hl, az ^ nht/hl, (3r ^ ht / {2hr) , Pz^ht/{2hz). (A3) 

For r — we have dT/dr = 0, therefore for small r one has T — a + O(r^), i.e. d^T/dr^ = dT/ [rdr). Hence, instead 
of (|A2p we have the following formula: 



To^i = a.(To"^-+i + + ^arT^^ + (1 - 2a. - Aaz)T^^ - f3zVzOjiT^,,+, - (A4) 

For z — we have T^g = Tq. 



3. Boundary interpolation for temperatures inside the drop 

Let G — dT/dn ~ —LJ{r)/k, and D{i,j) is the mesh point which is close to the surface. The point D is inside the 
drop and at least one of its nearest neighbors has to be outside the drop. In linear approximation T — a + br + cz 
in a vicinity of the point D(i,j). We denote the temperatures at points B{i — and C{i,j — 1) as ub and uc 
correspondingly. Then 

6 sin + c cos (/> = G, a + b{i — l)hr + cjh^ = ub, a + bihr + c{j — l)hz = uc- (A5) 

The solution of the set of equations is 

a = {hr coa(l){iuB — {i — l)uc) + hz aiii(/){uB + j{uc — Ub)) — hrhzG{j + i ~ 1)) / R (A6) 

b = {hzG + [uc - Ub) COB (i))/R (A7) 

c = {hrG + {uB - uc)s\n(t))/R (A8) 

R = hr cos (p + hz sin (p. (A9) 

The coefficients a, b, c allow to calculate the temperature at the point D{i, j) using the formula Tij = a + bh^i + chzj- 
Also, the coefficients a, 6, c allow to find the temperature at the surface of the drop near the point D. 



4. Alternating direction implicit method for vapor density 

In order to cover quite a large region for the vapor density calculation with Eq. (fT4|) . this is convenient to use 
irregular mesh. We use the following mesh: = ihr for i < rir, ri — hr{nr + 1.04'""''+^^ — 1.04^^) for i > Ur, Zj = jhz 
for j < Uz, Zj = hz{nz + 1.04J-"-+'^9 - 1.04^^) for j > n^. Here rir = Uz = 700, rirhr = R, Uzhz = h. This is the 
mesh with sufficiently small steps near the drop surface and with exponentially increasing steps, which are chosen in 
accordance with the asymptotic decay of the vapor density far away from the drop. 

We denote the distances between the mesh point and its nearest neighbors as a = — r^^i, b = ri+i — r^, 

c ~ Zj^i — Zj and d — Zj — Zj^i. Then the finite difference representations for the second derivatives are 
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= - (« + + (AlO) 

= ^^37-j('^<7+i - (c + ^X, + c<,_i), (All) 



where m^- is the vapor density at mesh point (i,j) for the n-th time step. 

We apply the ahernating direction method to Eq. (fT4| with the above notations. In the first part of the method 
one takes r-derivative imphcitly. Then the finite difference representation of Eq. p4p is 

n+l/2 „ n+1/2 n+1/2 

Z?V2 " +^^''v + r a + b ' ^"^^^^ 

For given vapor density at time step n it is convenient to rewrite this expression as 

+ « - 2/iDh,))ul+'/' + 4u^^f = c>^^_i + (rf; - 2/iDh,))u^^ + e>J^^+i. (A13) 

Here 

c', = 2/((a + b)a) - l/((a + &)r), = -2/(a6), e', = 2/((a + 6)6) + l/((a + 6)r), (AM) 

c;' = -2/((c + d)d), d;' = 2/(cd), e;' = -2/((c+d)c), (A15) 

c;, = 0, d^ = -4/a2, eo = 4/a^ (A16) 

c[,' = 0, d^,' = 2/c2, e[,' = -2/c^ (A17) 

For each j the tridiagonal matrix algorithm is used to solve the set of equations (|A13P for the vapor density at the 
time step n + 1/2. 

In the second part of the method one takes z-derivative implicitly and represents Eq. p4)) as 



n+l _ n+1/2 n+1/2 _ n+1/2 

Dht/2 -^^""^^ +r a + 6 ' ^^^^^ 

For given vapor density at time step n + 1/2 this expression takes the form 

c>i:+'i + (4 + ViDh,))ul+' + e';u^+l, = c^u:!^f + (d^ + 2/{Dh,))u-;'/' + ^^'^ , (A19) 

where the coefficients are given in (jAl4P - (|A17p . For each i the tridiagonal matrix algorithm is used to solve the set 
of equations (|A19p for the vapor density at time step n + l. 



5. The boundary interpolation for vapor density 

Consider a mesh point D{i,j) close to the surface. The point D is outside the drop and at least one of its nearest 
neighbors is inside the drop. In linear approximation u(r, z) = a + 6r + cz in a vicinity of the point D{i,j). Consider 



mesh points B{i + l,j), C{i,j + 1) and a point A of the drop surface near D. Then one has 

a + brA + czA — u{A), a + brs + czb — u{B), a + brc + czc — u{C). (A20) 

The solution of the set of equations is 

a = {{rcZB - rBZc)u{A) + (r^zc - rcZA)u[B) + {rsZA - rAZB)u{C))/ R, (A21) 

6 = {zc{u{A)-u{B)) + zb{u{C)-u{A)) + za{u{B)-u{C)))/R, (A22) 

c = {{tb - rc)u{A) + {rc - taMB) + - rB)u{C))/R, (A23) 

R = [rcZB ~ rBZc) + (rAZc - rcZA) + {tbza- tazb)- (A24) 



In the first part of the alternating direction method the calculations of rows proceed towards smaller values of j. For 
this reason one should consider here u{C) and u{A) ~ Ug as given quantities, whereas u{B) and u{D) are unknown. 
It is convenient under these conditions to represent (jA2ip - (jA24p as 



a = flo + aiu{B), b = bo + biu{B), c = cq + ciu{B), 



(A25) 
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and obtain explicit expressions for oq, ai, 5o, 61, cq, Ci. This results in linear relation between u{D) and u(B) 

u{D) = a + brn + czd = (ao + boro + co^d) + (oi + biro + cizi3)u(i?), (A26) 

which can be transformed to the form c'^u^^^^j^ +d'^u^j^^^^ + ~ b\. This completes the set of equations ()A13p 

for the tridiagonal matrix algorithm. Here c'^ = 0, = 1, e\ — —{a\ + h\rD + ciz/j), h\ = ao + b^VB + co^d. 

To carry out the boundary interpolation in the second part of the alternating direction method, similar expressions 
can be derived to relate m(-D) and u(C). 

The obtained values h and c also allow to find the local evaporation rate at the surface point A: 



J(A) = B 



= -£>(& sin (/) + c cos (/)). (A27) 



dn 

The total drop evaporation rate is found from local values (|A27p with numerical integration of . 

6. Alternating direction implicit method for stream function 

The equation to be numerically solved is (compare with p5p ) 

w at or^ oz^ r or 
The finite difference representation for the second derivatives on a regular mesh are 

5lipi;j = {'tjJi+i.j - 2i>ij + ipi^i,j)/hl, Sli>ij = - S'i/'y + tlJi,j-i)/hl. (A29) 

Taking r-derivative implicitly in the first part of the alternating direction method, we have for Eq. (jA28|) 

= , ^.,5 _ i*±M^^ - r-,„. ,A30) 

For a given ip at time step n it is convenient to rewrite this expression as 

c^V'I^f + « - 2/{i,htM;'/' + = + id- - 2/{u;hm^ + ,+1 + n^J. (A31) 

where 

= (1 + l/{2i))/hl, d[ = -2/hl e', = (1 - ll{2i))/hl (A32) 

c'; = ~l/hl d';^2/hl e';^-l/hl (A33) 

For each j the tridiagonal matrix algorithm is used to solve the set of equations (jA3ip and obtain ijj at time step 
n+ 1/2. 

In the second part of the method one takes 2;-derivative implicitly and represents Eq. (|A28|) as 
For given ■0 at time step n + 1/2 this expression takes the form 

c^V'- + K + 2/{u.h,m;^ + e';i^+i, = d^+^i; + « + 2i{uh,))i^i+"^ + ^^^^ti!" - ^In^ (ass) 

where the coefficients are defined in (jA32[) . (jA33[) . 

For each i the tridiagonal matrix algorithm is used to solve the set of equations (|A35P for the stream function ?/; at 
time step n + 1 . 
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7. Boundary interpolation for stream function 



Consider a mesh point D{i,j) close to the surface. The point D is inside the drop and at least one of its nearest 
neighbors is outside of the drop. In linear approximation ■(/'(r, z) — a + br + cz in a vicinity of the point D{i,j). For 
mesh points B{i — C{i,j — 1) and a point A of the drop surface near D one gets 

a + bvA + czA = 0, a + hvB + czb = ip{B), a + brc + czc = V'(C)- (A36) 

We solve the set of equations (jA36l) and obtain a, b, c as functions of ipiB), ip{C). The solution takes the form 
a = ao + aiipiyB), b = bo + bi4>iB), c = cq + ciip(B), where ao, bo, cq, ai, &i, ci are functions of ip{C). 

In the first part of the alternating direction method the calculations of rows inside the drop proceed towards larger 
values of J. For this reason one should consider here V'(C') ^-s a given quantity, whereas ipiB) and ipiD) are unknown. 
In order to complete the set of equations (jA3ip for the tridiagonal matrix algorithm we obtain the following relation 
between iJj{B) and ipiD) 

^Ij{D) = a + bra + czd = (ao + ^aro + cq-zd) + (ai + biro + ciZd)^{B). (A37) 

Similarly, for the calculation of columns in the second part of the alternating direction method, one obtains a = 
ao + aitp{C), b = bo + 6iV'(C), c = co + ci?/'(C) and 

■ipiD) = a + bro + czd (oq + 6o?'D + coZd) + (ai + biro + ciZd)'iI'{C). (A38) 



8. Explicit method for vorticity 



For solving Eq. ([5^ the following explicit finite difference method is used: 

In ''J I '^-^-1 I - f-:" _ 2-7" 4- -7" W 



^(7.':,+i - 27S- + Ti:,-!) + - 7r-i,) - ^7.,, (A39) 



+ (1 - 2a, - 2a, - ar/{2i^)h:, - PrV„j{l':+,,, - 7,"-!,^) " /3.«..j(7",+i 



7",-i)- (A40) 



Here 



vht/hl, a,^vht/hl, (3r ^ ht/ {2hr), (3, ^ ht/ {2h,). 



(A41) 



9. Boundary interpolation for vorticity 

Consider a mesh point D{i,j) close to the surface. The point D is inside the drop and at least one of its nearest 
neighbors has to be outside the drop. In linear approximation 7 = a + 6r + cz in a vicinity of the point D{i,j). We 
denote the vorticity at points B{i and C{i,j — 1) as 7(i?) and 7(C) correspondingly. The vorticity at the point 

A on the drop surface near the point D is obtained as 7(A) = da/{7]ds) + 2vrd(j)/ds (see Appendix [B|1 . Then 



a + brA + czA = 7(A), (A42) 

a + b{i-l)hr + cjh, = -f{B), (A43) 

a + bihr + c{j - l)h;, = -f{C). (A44) 

The solution of the set of equations is 

a = (h,rA{l{B) + (7(C) - 7(B))j) + Kza{i{C) + {j{B) - j{C))t) - j{A)hrK{t + j - 1))/R, (A45) 

b = (zAiiiC) - 7(i?)) + h^iliA) + 7(B)(j - 1) - niC)))/R, (A46) 

c = {rAh{B) - 7(C)) + hr{j{A) - tj{B) + (z - l)7(C)))/i?, (A47) 

R = h^rA + hrZA- hrh,{i + j - 1). (A48) 



The coefficients a, b, c allow to calculate 7 at the point D{i,j) using the formula j{D) — a + bhri + chzj. 
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10. Calculation of the drop surface 

For solving numerically Eqs.®, (fTO]) for the drop surface in the form of y(s) — (r(s), ^(s), z(s))-^, it is convenient 
to modify the initial conditions pT|) as follows: 

y(0) - (i?o<5,<5,0f . (A49) 

Here r{s), (j){s), z{s) are unknown functions, s e {0,Smax}, and S — 10~^ is introduced to fulfill the condition 
m/riO) = S/iRoS) = 1/Ro. 

Consider the set of points Sk = hk, where k = 0,1, . . . , K and Kh = Smax- Eqs. ([9]), (fTO|) . (|A49p are solved with 
the Runge-Kutta method: 

yfc+i = yfc + 5(pi+2p2 + 2p3 + p4), (A50) 
Pi - f(sfe,yfe), (A51) 

P2 = f (^Sfe + ^,yfc + ^Pi) , (A52) 

P3 = f (^Sfc + ^,yfc + ^P2^ , (A53) 
P4 = f (sfc + /i,yfc + /ip3). (A54) 

To increase the accuracy further, an interpolation is used for obtaining the values Smax and h. Then, the drop mass 
is calculated numerically with Eq. (|12p . and the evaporation rate is given by p^ . 

APPENDIX B: BOUNDARY CONDITIONS FOR VORTICITY AND STREAM FUNCTION 

For a derivation of the boundary condition at the surface of the drop for the quantity 7, which satisfies Eq. (I32p . 
we consider a vicinity of a point taken at the surface of the drop. The r- and z- components of the velocity can be 
expressed at the surface via the respective tangential and normal components, lying in the rz-plane: 

Vr = Wr COS (/) + U„ siu (j), (Bl) 

Vz — — Wr sin (/) + u„ COS 0. (B2) 

The angle 4> between r- and r- projections of the velocity depends, in general, on the coordinate s along the surface. 
Therefore, 

— = — — cos (p - Vrsmcj) ■ — + — — sm 0, (B3) 
oz oz az oz 

-TT— = — sm0 - t;^ COS0 • — + — — COS0. (B4) 
or or ar or 



It follows from Eqs. dBlj), 1^ and fBi 

dVr dVz ( dVr dVr .A f d(f> d(/) . \ ( dVn dVn 



7 = — — = — — cos (p + — - sm (p \ +Vr\ — COS - — sm - — — cos (p - — — sm = 

oz Or \ oz Or J \dr dz J \ Or Oz 

- ^ J^- y ^ - ^ - ^ + v — (B5) 

dn ds ds dn ds 
The boundary condition at the surface of the drop can be expressed as (45| 

where the unit vector n is directed towards the vapor along the normal to the surface. Taking the tangential component 
of this equation, we find 

da ( dvi dvk\ { dvi dvk 

ds \ dxk dxi J * \dn ^ ds 

- ri( — + - Vk^^^ - V ( ~v — ] (B7) 
\ dn ds ds J \ dn ^ ds ' 
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Here are the components of the unit vector r tangential to the surface. 

As follows from Eqs. (|B5p and (|B7[) . the boundary condition for Eq. ([5^ . i. e. for the quantity 7(7", z) at the surface 
of the drop, takes the following form 

1 da ^ dd> 

j=--^ + 2vr-^. B8 
7] ds ds 

Consider now the boundary conditions for ip, which satisfies Eq. (j35p . For z = we have i9?/'(r, z — 0)/dr = ~rvz{r, z — 
0) = and, hence, ip{r,z = 0) = const. At the symmetry axis r = we have dtpir = 0,z)/dz = rvr{r = 0, z) = 0, 
therefore ijj{r = 0, z) = const. At last, on the outer surface of the drop we have 

— — = — — cos — sm = —rVz cos — rvr sm = — rw„ = U. (B9j 

OS or Oz 

Therefore, ip — const on the surface of the drop, and, hence, ip — const at all boundaries. Since only the derivatives 
of the stream function tp enter expressions for physical quantities, the particular value of the constant does not have 
any observable consequences. For this reason, one can put "0 = throughout the boundary. 
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